% This code can be executed to generate Figure S4C. The blue-green curves
% are for the G1 pre-buckle (GG4; R=4 inch). The pink curves are for the G2
% pre-buckle (GG2; R=2.6 inch). The dashed lines show the pre-buckle factor 
% for the segment AE of bottom tape (GG2_ae and GG4_ae) which does not 
% change as a function of the junction angle (see Figure S4A and S4B). The
% solid curves are for the pre-buckle factor (G) of AP.
% Additional formatting is done externally in Inkscape.
% The folder titled, "Raw Data Test Files," must be open to run this code!

clear;
w=12.7;
t=0.058;
Rs=[66.04, 101.6];
% Initialize variables for complete computation
angles=linspace(20,90,1000);
APs=zeros(1,length(angles));
L2b=zeros(1,length(angles));
L2t=zeros(1,length(angles));
L4b=zeros(1,length(angles));
L4t=zeros(1,length(angles));
GG4=zeros(1,length(angles));
GG2=zeros(1,length(angles));
GG4_ae=zeros(1,length(angles));
GG2_ae=zeros(1,length(angles));
c=zeros(1,length(angles));

for j=1:length(Rs)
    R=Rs(j)+t/2;
    Rt=Rs(j)+1.5*t;
    for i=1:length(angles)
        %Initialize variables for particular junction angle
        phi=angles(i);
        c(i)=2.*R.*sin(w./(2.*R)).*cotd(phi);  
        
        omega=w./R;
        % parameters for bottom tape
        x=@(t) R.*cos(omega.*t);
        y=@(t) R.*sin(omega.*t);
        z=@(t) 2.*R.*sin(w./(2.*R)).*cotd(phi).*t;
        % parameters for top tape
        xt=@(t) Rt.*cos(omega.*t);
        yt=@(t) Rt.*sin(omega.*t);
        zt=@(t) 2.*R.*sin(w./(2.*R)).*cotd(phi).*t;
        APs(i)=w.*cscd(phi);
        % arc lengths of top (Lt) and bottom (Lb) tapes
        Lb=sqrt((-omega.*R).^2+(c(i)).^2);
        Lt=sqrt((-omega.*Rt).^2+(c(i)).^2);
        
        color=linspace(100,255,length(angles));
        if j==1
            L2b(i)=Lb;
            L2t(i)=Lt;
            %Determine pre-buckle factor G for segment AP
            GG2(i)=(Lt-Lb)./APs(i);
            %Determine pre-buckle factor for segment AE
            GG2_ae(i)=Rt./R-1;
            
            % Optional: Plot the actual curvature of segment AP
            %fplot3(x,y,z,[0,1],'Color',[color(i)/255, 0/255, 0/255]);
            %hold on
            %fplot3(xt,yt,zt,[0,1],'Color',[color(i)/255, 0/255, 0/255]);
        else
            L4b(i)=Lb;
            L4t(i)=Lt;
            %Determine pre-buckle factor G for segment AP
            GG4(i)=(Lt-Lb)./APs(i);
            %Determine pre-buckle factor for segment AE
            GG4_ae(i)=Rt./R-1;
        end
    end
end

%plot G as a function of junction angle for both radii
figure;
hold on
plot(angles,GG4,'Color',[0/255, 150/255, 150/255])
plot(angles,GG2,'Color',[255/255, 150/255, 150/255])
plot(angles,GG4_ae,'LineStyle','--','Color',[0/255, 150/255, 150/255])
plot(angles,GG2_ae,'LineStyle','--','Color',[255/255, 150/255, 150/255])

%Plot formatting
title('"Preload" Factors');
xlabel('Junction Angle \phi (\circ)','FontSize',14);
ylabel('G','FontSize',14);